#include "PyrolysisADMaterial.h"
registerMooseObject("TrilobitaApp", PyrolysisADMaterial);


InputParameters
PyrolysisADMaterial::validParams()
{
    InputParameters params = ADMaterial::validParams();
    
    params.addRequiredCoupledVar("temperature", "The temperature");
    params.addRequiredCoupledVar("porepressure", "The porepressure");
    return params;
}

PyrolysisADMaterial::PyrolysisADMaterial(const InputParameters& parameters)
    : ADMaterial(parameters),
    _temperature(adCoupledValue("temperature")),
    _temperature_old(coupledValueOld("temperature")),
    _porepressure(adCoupledValue("porepressure")),
    _grad_temperature(adCoupledGradient("temperature")),
    _grad_porepressure(adCoupledGradient("porepressure")),
    _rhos(getADMaterialProperty<Real>("solid_density")),
    _drhos_dt(getADMaterialProperty<Real>("drhos_dt")),
    _rhog(getADMaterialProperty<Real>("pyrolysis_gas_density")),
    _thermal_conductivity_s(getADMaterialProperty<Real>("solid_thermal_conductivity")),
    _specific_heat_s(getADMaterialProperty<Real>("solid_specific_heat")),
    _porosity(getADMaterialProperty<Real>("porosity")),
    _permeability(getADMaterialProperty<Real>("permeability")),
    _cpg(getADMaterialProperty<Real>("pyrolysis_gas_specific_heat")),
    _viscosity(getADMaterialProperty<Real>("pyrolysis_gas_viscosity")),
    _pyrolysis_heat(getADMaterialProperty<Real>("pyrolysis_heat")),
    _rhos_old(getMaterialPropertyOld<Real>("solid_density")),
    _rhog_old(getMaterialPropertyOld<Real>("pyrolysis_gas_density")),
    _thermal_conductivity_s_old(getMaterialPropertyOld<Real>("solid_thermal_conductivity")),
    _specific_heat_s_old(getMaterialPropertyOld<Real>("solid_specific_heat")),
    _porosity_old(getMaterialPropertyOld<Real>("porosity")),
    _cpg_old(getMaterialPropertyOld<Real>("pyrolysis_gas_specific_heat")),

    _Darcy_velocity(declareADProperty<RealVectorValue>("Darcy_velocity")),
    _mass_td_strong_residual(declareADProperty<Real>("mass_td_strong_residual")),
    _mass_viscous_strong_residual(declareADProperty<RealVectorValue>("mass_viscous_strong_residual")),
    _mass_source_strong_residual(declareADProperty<Real>("mass_source_strong_residual")),
    _energy_td_strong_residual(declareADProperty<Real>("energy_td_strong_residual")),
    _energy_td_strong_residual_factor(declareADProperty<Real>("energy_td_strong_residual_factor")),
    _energy_diffusion_strong_residual(declareADProperty<RealVectorValue>("energy_diffusion_strong_residual")),
    _energy_source_strong_residual(declareADProperty<Real>("energy_source_strong_residual")),
    _energy_convection_strong_residual_factor(declareADProperty<RealVectorValue>("energy_convection_strong_residual_factor")),
    _energy_convection_as_diffusion_strong_residual(declareADProperty<RealVectorValue>("energy_convection_as_diffusion_strong_residual")),
    _energy_convection_strong_residual(declareADProperty<Real>("energy_convection_strong_residual")),
    _use_displaced_mesh(getParam<bool>("use_displaced_mesh")),
    _ad_q_point(_bnd ? _assembly.adQPointsFace() : _assembly.adQPoints())
{
   
}

void
PyrolysisADMaterial::initialSetup()
{
   
}

void
PyrolysisADMaterial::computeQpProperties()
{
    ADReal mass = _porosity[_qp] * _rhog[_qp];
    ADReal mass_old = _porosity_old[_qp] * _rhog_old[_qp];
    _mass_td_strong_residual[_qp] = (mass - mass_old) / _dt;
    if ((_porosity[_qp]* _viscosity[_qp]) > 0)
    {
        _Darcy_velocity[_qp] = -_permeability[_qp] * _grad_porepressure[_qp] / _viscosity[_qp];
    }
    else
    {
        _Darcy_velocity[_qp] = 0 *_grad_porepressure[_qp];
    }
    _mass_viscous_strong_residual[_qp] = -_rhog[_qp] *  _Darcy_velocity[_qp];
    _mass_source_strong_residual[_qp] = _drhos_dt[_qp];
    ADReal energy_solid = _rhos[_qp] * _specific_heat_s[_qp] * _temperature[_qp];
    ADReal energy_gas = _porosity[_qp] * _rhog[_qp] * _cpg[_qp] * _temperature[_qp];
    ADReal energy_solid_old = _rhos_old[_qp] * _specific_heat_s_old[_qp] * _temperature_old[_qp];
    ADReal energy_gas_old = _porosity_old[_qp] * _rhog_old[_qp] * _cpg_old[_qp] * _temperature_old[_qp];
    _energy_td_strong_residual[_qp] = ((energy_solid+ energy_gas)-(energy_solid_old+ energy_gas_old)) / _dt;
    _energy_td_strong_residual_factor[_qp] = (_rhos[_qp] * _specific_heat_s[_qp] + _porosity[_qp] * _rhog[_qp] * _cpg[_qp]);
    _energy_diffusion_strong_residual[_qp] = _thermal_conductivity_s[_qp] * _grad_temperature[_qp];
    _energy_source_strong_residual[_qp] = -_pyrolysis_heat[_qp]*_drhos_dt[_qp];
    _energy_convection_strong_residual_factor[_qp] = _porosity[_qp] * _rhog[_qp] * _cpg[_qp] * _Darcy_velocity[_qp];
    _energy_convection_as_diffusion_strong_residual[_qp] = _porosity[_qp] * _rhog[_qp] * _cpg[_qp] * _temperature[_qp] * _Darcy_velocity[_qp];
    _energy_convection_strong_residual[_qp] = _porosity[_qp] * _rhog[_qp] * _cpg[_qp] * _Darcy_velocity[_qp] * _grad_temperature[_qp];
   // std::cout << _energy_convection_strong_residual[_qp].value() << std::endl;
}


